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ABSTRACT 

We present results from a suite of axisymmetric, core-collapse supernova simulations 
in which hydrodynamic recoil from an asymmetric explosion produces large proto- 
neutron star (PNS) velocities. We use the adaptive-mesh refinement code CASTRO 
to self-consistently follow core-collapse, the formation of the PNS and its subsequent 
acceleration. We obtain recoil velocities of up to 620 km at ~1 s after bounce. These 
velocities are consistent with the observed distribution of pulsar kicks and with PNS 
velocities obtained in other theoretical calculations. Our PNSs are still accelerating 
at several hundred kms~^ at the end of our calculations, suggesting that even the 
highest velocity pulsars may be explained by hydrodynamic recoil in generic, core- 
collapse supernovae. 

Key words: supernovae: general - hydrodynamics - stars: interiors - pulsars - 
neutron stars 



1 INTRODUCTION 

At birth, pulsars achieve velocities well above those of their 
progenitor population ( Gunn and Ostriker||1970 Lyne and 



Lorimer 



19941 

T 



These pulsar "kicks" typically range from 
200 kms~^ to 400 kms~^, with the fastest neutron stars 
achieving velocities near, or in excess of, 1000 kms~^ (|Lyne| 



and Lorimer|1 994':'Chatterjee et al.'2 005[[Hobbs et al.|2005 l 
The current distribution of observed pulsar velocities is 
Maxwellian, hinting at a common acceleration mechanism 



2005 Zou et al. 



(Hansen and Phinney 1997 Hobbs et al 
|2005[ [Faucher-Giguere and Kaspi|2006| ). 

Many scenarios have been proposed for the origin of 
pulsar recoil and neutron star kicks. Popular mechanisms 
often require strongly magnetized systems, exotic neutrino 
physics, and/or rapid rotation to produce substantial kicks. 
For example, in the presence of strong magnetic fields, 
neutrino-matter interactions can generate neutron star ve- 
locities on the order of a few hundred kms~^ by producing 
1% dipole asymmetries ( Lai and Qian|1998 Nardi and Zu- 



luaga2001"Lai et al.'200lVKuse nko and Se gre'1999:'Lambl 
ase.2005j Barkovich et al. 2004; Fuller et al.,,2003 ; Kishimoto 



2011 ). Many of these scenarios require magnetic fields in the 
magnetar range (i.e. IQ^*"^^ G) and may not produce sub- 
stantial kicks in typical core-collapse supernovae. Other sce- 
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narios involve jet/counter-jet misalignment launched near 
the proto-neutron star (PNS). In such situations, the jets 
accompany an associated gamma-ray burst (GRB) or form 
through magneto-rotational processes during core collapse 



Ccn 1998; Khokhlov et al. ,1999; Sawai et al. 2008; Pa^ 



tCc 

psh and Soker||2011[ ). These scenarios require rapid rota- 
tion and therefore, may only manifest in a small subset of 
core-collapse events. 

If neutron star kicks are a generic outcome of core col- 
lapse, then a natural explanation is recoil during an asym- 
metric supernova explosion. In the current, most sophis- 
ticated simulations, the bounce shock, launched when the 
equation of state stiffens at nuclear densities, stalls due to 
thermal energy losses from neutrino emission and dissocia- 
tion of nuclei into nucleons. The stalled shock itself is subject 
to hydrodynamic and neutrino-driven instabilities, which 
manifest as prominent low-order ^-mode oscillations in ax- 



isymmetric simulations of non-rotating progenitors ( Blondin 
et al.||2003| [Scheck et al.|[200^ 



et al. 2006 , Burrows et al. |2007a 



Buras et al. 



2006 



Scheck 



JBlondii^^anc Mczzacappa 



2007j|Ott et al.|2008nFerndndez and Thompso'nr2009t,Nord' 



haus et al.|2010b|a||Fernandez|2010[ [Brandt et al.|2011[ ). At 

the onset of shock revival, the PNS may recoil if large-scale 
asymmetries are present during the ensuing supernova ex- 
plosion. While the mechanism by which core-collapse super- 
nova progenitors explode is not fully understood, the most 
probable scenario involves absorption of neutrinos in the 



post-shock "gain region" ( Bethe and Wilson 1985 1 and likely 
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requires the development of multi-dimensional instabilities 
in fully three-dimensional radiation-hydrodynamic simula- 
tions to succeed (Nordhaus et al. 2010b I. Nonetheless, re- 



coil from an asymmetric neutrino-driven explosion presents 
a natural mechanism by which PNSs achieve high velocities 



( Burrows and Hayes||1996| [Scheck et al.p004 |2006| [Nord- 
haus et al.|2010a Wongwathanarat et al.|2010[ ) and appears 
to be supported by recent X-ray observations of the Cas- 
siopeia A supernova remnant ( Hwang and Laming|20lT l. 

Computational studies of recoil require multi- 
dimensional, radiation-hydrodynamics calculations which 
start at the onset of collapse and follow the dynamics 
self-consistently. This includes the formation of the PNS, 
the explosion, the propagation of the shock front through 
the stellar envelope and eventually, decoupling of the 
PNS from the surrounding material. Such an approach is 
computationally challenging and as such, various techniques 
have been adopted to make the computations tractable. 
One popular approach is to commence the calculations 
after bounce by mapping spherically symmetric solutions 
onto a multi-dimensional grid and excising the PNS from 
the computational domain (Scheck et al. 2004 2006 



Wongwathanarat et al. 20101. This requires one to infer a 



PNS kick through a rigid, impenetrable inner boundary, 
but allows one to track the supernova explosion for several 
seconds and to distances greater than 10,000 km. While this 
approach is appealing, it must be checked by simulations 
which include the PNS in the computational domain. 

Recently, we have carried out the first axisymmetric, 
radiation-hydrodynamic simulation of recoil with the multi- 
group, arbitrary, Lagrangian-Eulerian code VULCAN/2D 
( [Nordhaus et al.||2010a| . By transitioning from a spherical- 
polar mesh to a pseudo- Cartesian mesh at the center of the 
domain, we self-consistently tracked the PNS's formation 
and off-center motion. This calculation was computationally 
expensive, as it employed multi-group flux limited diffusion 
neutrino transport and followed the supernova explosion un- 
til the shock reached the 5,000 km radial outer boundary 
of the domain. At that time, the PNS had reached a ve- 
locity of ~150 kms~^ but had yet to fully decouple from 
the ejecta. The PNS recoil was due almost entirely to hy- 
drodynamical processes and was consistent with previous 
excised-core calculations ( Scheck et al.[[2004 2006 1 . Asym- 



metric neutrino emission contributed ~2% to the overall kick 
magnitude. This suggests that neutrinos play no significant 
role in accelerating neutron stars to high velocities during 
typical core-collapse supernovae. At the end of the calcula- 
tion, significant acceleration (~350 kms~'^), in addition to 
the degree of asymmetry in the ejecta, suggested that higher 
PNS velocities were possible. Verifying these estimates re- 
quires tracking the supernova shock to greater radial dis- 
tances and later times. 



To expand upon the work of Nordhaus et al. (2010a I, 
we carry out a suite of axisymmetric collapse calcula- 
tions with the adaptive- mesh-refinement (AMR), radiation- 
hydrodynamics code, CASTRO ( [Almgren et al.|2010|[Zhang 
et al. 20111. By employing AMR and a simplified form of 



transport (see Sec.[2|, we expand the domain to a radial dis- 
tance of 10,000 km and perform multiple calculations. Our 
PNSs achieve recoil velocities ranging from tens of kms~^ 
up to ~620 kms~^. In general, the magnitude of the recoil 
depends on the degree of asymmetry at the time of explosion 



and the energy of the explosion itself ( Burrows et al.|2007b I . 
In Sec. 3, we discuss the physical processes that accelerate 
the PNS. In Sec. |4j we compare our results with previous 
work before concluding and discussing future work in Sec. 

m 



2 COMPUTATIONAL METHODOLOGY 

We carry out our simulations using the AMR, radiation- 
hydrodynamics code, CASTRO ( [Almgren et al.|2010|[Zhang 
[et al.[[201l| ). CASTRO is a finite- volume, Eulerian code 
which employs an unsplit version of the piecewise parabolic 
method (PPM) for the hydrodynamics and a multigrid Pois- 
son solver to handle self-gravity. To facilitate multiple calcu- 
lations, we simplify our radiation transport by using radius- 
and temperature-dependent prescriptions for the neutrino 
heating and cooling rates H and C. We solve the fully com- 
pressible equations of hydrodynamics: 



dtp = -V ■ (pv) 
dt (pv) = -V ■ (pvv) - Vp + pg 
dt (pE) = -V ■ {pvE + pv) + pv ■ g + p{H - C) , 



(1) 



where p, T, p, g, and v are the fluid pressure, temperature, 
density, gravitational acceleration, and velocity. The specific 
total energy is given as E = e + ^v^ , where e represents the 
internal energy. Neutrino heating and cooling occur via the 
super-allowed charged-current reactions involving free nucle- 
ons, electrons, protons, electron neutrinos and anti-electron 
neutrinos. We use the heating and cooling rates derived in 
[Janka[ ( [2001|l and previously used by [Murphy and Burrows[ 



(2008) and [Nordhaus et al.[ ([2010b[). These rates, assuming 
the electron and anti-electron neutrino luminosities, L^^ and 
Lj}^, to be equal, are 



U = 1.544 X 10^ 



1052 erg s-iy V4 MeV 



(2) 



^ 100km y 



erg 



g s 



and 



C = 1.399 X 10^ 



2 MeV 



(K + Vp) e 



erg 



g s 



(3) 



where T^^ is the electron neutrino temperature, r is the dis- 
tance from the center of the star, Yn and Y-p are the neutron 
and proton fractions, and Teff is an effective optical depth for 
electron and anti-electron neutrinos (see Eqs. 6 & 7 of[Hanke[ 



et al. 



(20111). The factor e effects a transition between 



the dense inner regions, where neutrinos are trapped and 
heating and cooling are suppressed, to the outer, optically 
thin regions. To avoid the prohibitive cost of global calcula- 
tions, we fit Teff as a function of density by post-processing 
models and applying Eq. 6 of [Hanke et al.[ ( [2011| |. In the 
transition region, 10~^ < Tcfi < 1, the resulting power-law 
fits give root-mean-square residuals of ~10% in e~'^''^^ . Note 
also that the bulk of neutrino heating occurs within the in- 
ner few hundred kilometers and falls off rapidly with radial 
distance. The electron fraction, Y^, evolves on infall via a 
density prescription give n in [Liebendorfer et al.[ ([2005[ ) and 
previously employed by Murphy and Burrows ( 2008 j and 
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Nordhaus et al. (2010b I. To close the equations, we use the 



sophisticated, finite temperature, nuclear equation of state 



of Shen et al. ( 1998 1 



The current theoretical consensus holds that most 2D 
core-collapse simulations do not produce neutrino-induced 
explosions for the majority of progenitors without supple- 
menting the neutrino luminosity ( Burrows et al.|2007b| |Ott 



et al.|200 8"Nord haus eraL| r2010a: B randt et al.|2011||Fuj 



moto et al.^2011) . However, explosions in axisymmetric sim- 
ulations have been obtained after ~600 ms for a 15 Mq pro- 
genitor when a soft nuclear equation of state and variable 



Eddington factor closure technique are employed (Marek 
[and Janka||2009] ). We therefore induce explosions by vary- 
ing the neutrino driving luminosity Lv, (= La^) from 2.1 to 
2.5 X 10^^ ergs~^. These values are comparable to the time- 
dependent values achieved in 2D calculations with detailed 
neutrino transport ( Ott et al.|2008 Brandt et al.|20lT K Our 
simplified transport scheme couples this energy somewhat 
more efficiently to the post-shock region, and our results 
should be checked by future studies that employ sophisti- 
cated, but computationally expensive, neutrino transport. 
We hold the driving luminosity L^^ fixed in each calcula- 
tion, but vary it from run to run. Our simplified transport 
method, in conjunction with AMR, allows us to follow col- 
lapse, formation of the PNS, and any subsequent accelera- 
tion due to the supernova explosion. 

We use the 15- Mq, red-supergiant, non-rotating, solar- 



metallicity progenitor of Woosley and Weaver (19951. Our 



simulations begin at the onset of core-collapse, continue 
through the formation of the PNS and the launching, 
stalling, and revival of the bounce shock, and end when the 
shock approaches the edge of the computational domain. We 
employ a 2D, 10,000-km by 20,000-km domain discretized 
with a uniform coarse grid of 64 by 128 cells covering the 
fufi 180° range in polar angle. Within 200 km of the PNS, 
we use 4 AMR levels, each increasing the resolution by a fac- 
tor of 4, giving a minimum cell size of ~0.3 km. Exterior to 
200 km, the adaptive mesh tracks the high entropy, shocked 
material. 

We index our eight simulations by their driving electron 
neutrino luminosity L^^ in units of lO^'^ ergs~^, from L2.1 
to L2.5. For each simulation, we follow the explosion for ~1 
second of post-bounce evolution. For all but the L2.1 model, 
the calculation ends when the shock approaches the edge 
of the computational domain, 10,000 km from the PNS. At 
this point the total momentum on the grid is conserved to 
within ~1% of the core's final value. 



3 ACCELERATION OF THE PNS 

The hydrodynamic fiow behind the stalled shock is tur- 
bulent, and is soon deformed by the development of low- 



mode hydrodynamic and neutrino-driven instabilities ( Bur-| 



rows et al.|1995i|Blondin et al.|2003||Scheck et al.|2004|[2006| 
Blondin and Mezzacappa|2007 Fernandez|2010 1 . This guar- 
antees an asymmetric distribution of material when neutrino 
heating revives the stalled shock. The ensuing asymmet- 
ric explosion accelerates the PNS as the shock propagates 
through the stellar envelope. 

To understand the physical processes governing the re- 
coil, we first present the PNS velocities obtained in our sim- 
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Figure 1. PNS recoil velocities as a function of time after forma- 
tion of the bounce shock. The subscript in the simulation names 
refers to the electron neutrino, Ue, and anti-electron neutrino, Pe, 
luminosities, each in units of lO^'^ erg (see Eqs. ([2J1 and ([sjl 
and Table [1}. 



ulations. We compute the core positions as the centroids 
of the density distributions and differentiate to obtain the 
velocities, which we show as functions of time in Figure [l] 
The subscript in each simulation label indicates the driving 
L^^ in units of 10^^ ergs~^. The value of the driving Lj?^ 
is taken to be the same, thus yielding a total driving lumi- 
nosity in each simulation of itot = L,^^ + L^^ — 2L^^ . We 
obtain PNS velocities ranging from tens to many hundreds of 
kms~^. The highest velocity PNS (solid red curve in Fig.[l]) 
reaches a speed of 624 km s~^ after ~800 ms of post-bounce 
evolution, and is still accelerating at ~1000 kms~^. 

Figure [2] depicts the evolution of a representative run 
(model L2.4). The left column presents the evolution of elec- 
tron fraction, Yc, over the inner ~200 km, as the explo- 
sion progresses. The middle of the computational domain is 
marked by a black line {Z = 0). The PNS is clearly visi- 
ble as the deleptonized blue region, which begins a.t Z — 
and moves off-center. The right column shows the global 
evolution of the anisotropic supernova shock, with the color 
map depicting entropy. The explosion occurs primarily in 
the ~Z direction while the PNS recoils in the +Z direction. 
By ~800 ms after bounce, the PNS has largely decoupled 
from its surroundings. The middle column shows the density 
distribution with the shock outlined in pink. A high-density 
region above the PNS combines with a low-density region 
below it to gravitationally accelerate the PNS in the +Z 
direction. 



3.1 Physics of the Recoil 

To determine the physical processes governing the move- 
ment of the PNS, we post-process our results by comput- 
ing the hydrodynamic acceleration Sc of the core due to 
anisotropic pressure forces, momentum flux, and gravita- 
tional forces. The Eulerian equations of hydrodynamics give 



Gr 



dm - 



1 



PdS + 



pv^vdS 



(4) 
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Figure 2. The left column shows the evolution of the inner ~200 km of the simulation domain for model L2.4. The color map depicts the 
electron fraction, with velocity vectors overlaid; the black lines indicate Z = 0. The color map for the middle column depicts density. 
The pink curve shows the location of the shock. Exterior to the shock, the flow is radially inward. The right column panels show the entropy 
evolution of the explosion. The steep entropy jump just interior to the shock depicts the region where nucleons are reassociated into nuclei 
and alpha particles. At 707 ms after bounce, the PNS wind is seen as the dark blue region interior to the shock in the middle panel. The 
supernova explosion primarily occurs in the —Z direction, while the PNS recoils in the +Z direction. At ~800 ms after bounce, the PNS 
has largely decoupled from the surrounding material, but is still being accelerated by the gravitational pull of slow-moving ejecta in the 
+Z direction (see Fig. [sj. 



where p is the density, Mc and Vc are the mass and mean 
velocity of the core, P is the gas pressure, v is the fluid 
velocity, v-r is the radial component of the velocity, and rc 
is a fiducial spherical radius ( [Scheck et al.|[2006, ^Nordhaus] 
|et al.|2010a"| . 

The three terms in Eq. Q represent the contribu- 
tions to the acceleration from the gravitational field of mat- 
ter exterior to rc, anisotropic gas pressure, and momen- 
tum flux through rc, respectively. The first term assumes 
a spherically-symmetric matter distribution interior to rc, 
which in practice is an excellent approximation. In the limit 
of an isotropic explosion, no acceleration occurs and each 
term vanishes individually. 
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Equation ([4| includes contributions from hydrodynamic 
processes, but neglects radiation pressure asymmetries, 
which are not captured by our heating and cooling prescrip- 
tion (Eqs. |[2|-(|3|). In our previous calculations, which per- 
formed radiative transfer using multi-group flux-limited dif- 
fusion, neutrino momentum contributed < 2% of the overall 



kick (Nordhaus et al. 2010a I. This is consistent with previous 



studies that found neutrino radiation pressure to contribute 
~5% to the final kick velocity ( |Scheck et al.|2004| [2006( 1 . 

The relative contributions of the terms in Eq. ([4|) de- 
pend on the properties of the flow and the explosion dynam- 
ics. As a consequence of the explosion, the PNS generically 
recoils away from the high- velocity ejecta and towards the 
lower- velocity ejecta. Ifowever, the interpretation of the kick 
is not as straightforwar d as Eq. ([4| woul d suggest. As dis- 
cussed in Section III of Nordhaus et al. (2010a I, pressure 



and gravity do work on fluid elements; anisotropic pressure 
or gravitational forces at a small value of will appear as 
anisotropic momentum flux at a larger value of r^- 

Using a fiducial radius of = 200 km and integrating 
Eq. [4]from core bounce. Figure [3] shows the decomposition 
of the PNS kick into three components for models L2.25 (top 
panel), L2.3 (middle panel) and L2.4 (bottom panel). Note 
that the velocities have been reflected for model L2.3. In each 
panel, the dash-dotted black curve represents the smoothed 
centroid velocity, while the solid red curve is the sum of the 
three terms in Eq. Q. The gravitational component (short- 
dashed green curve) dominates the late-time evolution in 
all three simulations. The anisotropic pressure term (dot- 
dashed pink curve) flattens towards the end of each run as 
the PNS decouples from the ejecta. 

Model L2.3 achieved the highest kick velocity in our 
suite of simulations, more than 620 kms"^ at ~800 ms after 
bounce. The middle panel of Figure [3] shows its velocity evo- 
lution, decomposed using Eq. Q and inverted to facilitate 
comparison with models L2.25 and L2.4. Anisotropic pres- 
sure and momentum flux (dot-dashed pink and long-dashed 
blue lines, respectively) contributed almost nothing to the 
kick after ~400 ms from core bounce. Driven by the gravi- 
tational term in Eq. (Ml , this model was still accelerating at 
more than 1000 km s~ when the shock reached the edge of 
the computational domain. 

While they did not achieve as large a PNS velocity as 
L2.3, models L2.25 and L2.4 were still accelerating at ~1000 
and ~600 kms~^, respectively, at the end of our calcula- 
tions. In both cases, and particularly for the L2.4 model, 
this acceleration was dominated by the gravitational term in 
Eq. (|4|. Figure [2] clearly shows the PNS and ejecta in model 
L2.4 decoupling at ~650 ms after bounce (second panel from 
bottom), and having almost completely decoupled by ~800 
ms after bounce. 

The three models presented here comprise a represen- 
tative sample of our simulation results. Table [l] presents 
additional information on each of our runs, including the 
velocity, the explosion energy iJexp, and a, a dimensionless 
measure of the degree of asymmetry, at the end of the calcu- 
lations. The parameters a and -Eoxp are defined by Eqs. (|5| 
and ([6| in the following section. All of our calculations, ex- 
cept for L2.1, ended when the shock approached a radius of 
10,000 km; model L2.1 ended with 7?ahock = 3300 km. For 
a detailed discussion of the limitations of our approach, the 
effect of fixing the neutrino luminosities and the reliability of 
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Figure 3. Decomposed PNS kick velocities for models L2.25 (top 
panel), L2.3 (middle panel) and L2.4, obtained by integrating 
Eq. from bounce. We have inverted the velocities for the L2.3 
run. The dash-dotted-black curve depicts the PNS velocity com- 
puted using the centroid of the density, while the solid red curve 
shows the contributions to the kick from momentum flux (long- 
dashed blue curve), gravity (dashed green curve) and pressure 
(dash-dot pink curve). In each run, the PNS is still accelerating 
at more than 500 kms~^ at the end of the calculation, and this 
acceleration is dominated by the gravitational term. 



the late-time acceleration and explosion energies see Section 
[331 



Asymmetries in the Ejecta and Explosion 
Energies 



3.2 



The acceleration of the PNS depends on the dynamics of the 
explosion and the evolution of the asymmetry of the shocked 
material. This asymmetry may be quantified in various ways 
( Scheck et al.|2006 Burrows et al.|[2007b l; here, we adopt 

«=(«.>/{|^l> , (5) 



where {) denotes a mass-weighted average over the post- 
shock region with r > 100 km (thereby excluding the PNS 
itself). This choice is similar to the a presented in Scheck 
et al.| ( |2006l). 

Figure |4| shows the evolution of a for models L2.4 (solid 
blue curve) and L2.5 (dashed blue curve). The solid red and 
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Table 1. Parameters at the end of the simulation, when -Rghock ~ 
10,000 km; a and -Eoxp are calculated using Eqs. ([sjl and 
respectively. Model L2.1 ended with iishock — 

3300 km. Note 

that the PNS wind contributes ~50% of the explosion energies 
listed below. 



Model 


■upNS [km s 1] 


E^^p [1051 erg] 


a 


L2.1 


-40 


0.29 


0.026 


L2.15 


212 


0.69 


-0.25 


L2.2 


-186 


0.89 


0.08 


L2.25 


315 


0.69 


-0.23 


L2.3 


-624 


1.13 


0.23 


L2.35 


194 


1.28 


-0.06 


L2.4 


431 


1.23 


-0.15 


L2.5 


276 


0.99 


-0.10 



dashed red curves depict the PNS recoil velocities for the 
L2.4 and L2.5 models respectively. Our suite of simulations 
produced final values of a between —0.25 and 0.25 (see Ta- 
ble [l}. Since momentum is conserved, larger values of a lead 
to larger PNS recoil velocities. 

Figure [6] shows the position of the shock at the end of 
the calculation for five of our models; models with a neg- 
ative kick have been refiected. While the shock asymmetry 
does correlate with the kick velocity, the magnitude of the 
kick depends on the distribution of matter behind the shock, 
which we paramtrize using a. 

The third column of Table [T] shows the explosion energy 
at the end of the simulation, defined as the total energy of 
all unbound material on the grid. 



/ 

JEt 



P Mint + 



G 



Mo, 



(6) 



'Btot>0 

where Mint is the specific internal energy and A/enc is the 
mass interior to the fiuid element. At the end of our calcu- 
lations, the internal energy of shocked material dominates 
the kinetic energy by a factor of ~3-5, and the explosion 
energy is still increasing due to sustained neutrino heating 
(Eq. ([2|). The bulk of this heating occurs primarily within 
the first few hundred kilometers of the PNS and is driving 
the late-time PNS wind. 

The internal energy of shocked material will ultimately 
be converted into kinetic energy by pdV work. In this limit, 
the PNS kick will be a function of the anisotropy of the 
ejecta and -Ecxp- While the explosion will be nearly spheri- 
cal in the outer envelope, the anisotropy in the inner mass 
shells should freeze out at values close to those indicated 
in Table 1. This anisotropy in the inner ejecta velocities, 
with the bulk of the ejecta traveling opposite to the recoil- 
ing PNS, should be observable in the supernova remnant 
and is a specific prediction of our model. 

Figure |5] shows the total explosion energy (solid blue 
curve) as a function of time for model L2.4. The explosion 
energy of material in the bottom half of the computational 
grid {Z < 0) is depicted by the dot-dashed green curve, 
while the explosion energy of material in the top half of the 
computational grid {Z > 0) is depicted by the dashed red 
curve. Consistent with Fig. [2] this demontrates that model 
L2.4 explodes primarily in the —Z direction. The bottom 
panel of Fig. [5] shows the PNS mass as a function of time 
(solid blue curve) and the total mass of shocked, bound ma- 
terial exterior to the core (dot-dashed red curve). By the 
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Figure 4. The blue curves show the evolution of alpha (Eq. |5|) 
for models L2.4 (solid) and L2.5 (dashed), while the red curves 
show the core recoil velocity as a function of time. The core veloc- 
ity is always opposite the ejecta asymmetry due to conservation 
of momentum. 
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Figure 5. The top panel shows the explosion energy (Eq. ^) as 
a function of time for model L2.4. The solid blue curve depicts 
the total explosion energy. The explosion energy of material with 
Z < is shown by the dot-dashed green curve while the explosion 
energy of material with Z > is shown by the dashed red curve. 
The bottom panel presents the PNS mass (solid blue curve) as 
a function of time; its final value is 1.45 Mq. The dot-dashed 
red curve in the bottom panel shows the total shocked, bound 
mass exterior to the core as a function of time. At the end of the 



simulation, only • 
remains bound. 



TO Mq of shocked material outside the PNS 



end of model L2.4, the PNS mass is 1.45 Mq while there is 
little shocked, bound matter outside the PNS itself. 



3.3 Late-Time Evolution 

Our simplified transport scheme allows us to perform multi- 
ple calculations of the inner regions of a core-collapse su- 
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Figure 6. The position of the outer edge of the gain region toward 
the end of the calculations. Runs with negative PNS velocity have 
been reflected in this plot. The arrows depict the core velocity for 
each run (right axis). 



pernova. To date, only one time-dependent, multi-group, 
flux-limited diffusion, neutrino transport kick calculation, 
which includes the PNS on the computational grid, exists 
in the literature ( Nordhaus et al.|2010a ). As computational 
methods and resources improve, it will become possible to 
self-consistently connect the PNS kick with the large-scale 
properties of the ejecta. In this section, we discuss the late- 
time evolution of our simulations and some of the limitations 
of our approach. 

As previously mentioned, the magnitude of the PNS 
kick will increase with the degree of asymmetry in the ejecta 
and the explosion energy of the supernova. At the end of our 
calculations, the shocked matter's internal energy exceeds its 
kinetic energy by a factor of ~3-5. Adiabatic expansion will 
convert this internal energy into kinetic energy as the shock 
propagates through the stellar envelope. 

Our constant driving L^^ also deposits energy into the 
expanding ejecta, both by neutrino absorption in the gain 
region and by driving a ~0.1 Mqs~^ wind from the sur- 
face of the PNS. This late-time wind contributes ~50% of 
the explosion energy at the end of each simulation and typi- 
cally contains ~0.05 Mq of material. While neutrino-driven 



winds from the PNS are expected (Burrows et al. 19951, 



future improvements to this work should include more so- 
phisticated transport approaches which naturally incorpo- 
rate time-variable neutrino luminosities. 

We tested the effect of decaying neutrino driving lumi- 
nosities by restarting model L2.3 400 ms after bounce with 
an exponentially decreasing driving luminosity. We used an 
exponential decay timescale of 1 second, giving a ~35% re- 
duction in Lv^ at the end of our calculation. As a result of 
the lower energy injection rate, the late-time PNS wind de- 
creased by nearly a factor of two, the PNS took longer to 
decouple from the post-shock material, and the final PNS 
velocity was ~25% lower. Still, the PNS was accelerating 
gravitationally at ~1000 kms~^, nearly as fast as in the 
model with a constant driving luminosity. 



We also note that our calculations end when the shock 
reaches a fixed radius of ~10,000 km, rather than after a 
fixed amount of post-bounce time. The total amount of en- 
ergy injected into our models thus varies widely, making it 
difficult to connect the derived explosion energies with our 
observed kick velocities. 



4 COMPARISON WITH PREVIOUS 
NUMERICAL WORK 

Previous computational studies of pulsar recoil have em- 
ployed various simplifications and approximations to make 
the calculations tractable. These approaches include excis- 
ing the PNS from the computational domain, starting the 
calculations ~20 ms after bounce, and employing simpli- 



fied neutrino transport schemes (Scheck et al. 2004 2006 



[Wongwathanarat et al.||2010| ). The exclusion of the PNS 
from the domain is particularly useful, as it avoids the severe 
Courant limitation imposed by resolving the PNS. The PNS 
is replaced by a rigid, spherical boundary, which contracts 
according to a prescription from a detailed spherically- 
symmetric-collapse calculation. This approach is designed 
to mimic the settling of material as the PNS cools. 

By using all three of these simplifications, previous stud- 
ies have been able to track the shock to large dis tances 
(> 10"^ km) and late times (>1 s) ([Scheck et al.j2004 



2006 



Wongwathanarat et al. 20101. While useful for calculating 



long-term evolution, this approach requires inferring recoil 
through a rigid boundary of infinite inertial mass. Further- 
more, this approach neglects effects resulting from the dis- 
placement of the PNS relative to the surrounding matter. 
To compensate for this effect, these authors have added a 
kick to the gas which mimics movement of the PNS. The 
physical fidelity of such approximations has been verified by 
self-consistent calculations such as those presented in ( |Nord-| 
haus et al.|2010a l and in this work. 



Recently, Nordhaus et al. (2010a I presented the first 



axisymmetric, multi-group, flux-limited diffusion neutrino 
transport calculation of recoil in which core collapse lead to 
significant acceleration of a fully- formed PNS. The authors 
used the multi-group, arbitrary Lagrangian-Eulerian (ALE), 
radiation-hydrodynamics code VULCAN/2D ( Livne|1993 1. 
The calculation employed multi-group flux-limited diffusion 
neutrino transport (Livne et al. 20041, supplemented the 



neutrino luminosity by an additional = = 2 x lO^'^ 
ergs^^, and used the same 15- Mq progenitor as this work. 
During the induced, neutrino-driven explosion, a ~10% 
anisotropy in the ejecta led to a PNS recoil velocity of ~150 
kms~^ at the end of the calculation, when the shock reached 
a radius of ~5,000 km. Such a result in terms of PNS ve- 
locity and ejecta asymmetry compares favorably with model 
L2.2 presented in this work (see Table [T]) and the results of 



Scheck et al. ( 2006 1 



In general, given the different computational techniques 
and the use of three different codes, the agreement between 
our results and those of lScheck et al.l and INordhaus et al.l 
( |2010a| ) is gratifying. Our detailed calculations of the first 
second of post-bounce evolution produce high-velocity re- 



coils comparable with those in Scheck et al. ( 2004 2006 1 and 
[Nordhaus et al.| ( |2010a| l while following the evolution of the 
PNS itself. Taken together, these studies strongly suggest 
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that asymmetric core-collapse supernovae naturally lead to 
acceleration of the PNS and are capable of birthing the high- 
est velocity pulsars. 

While axisymmery would restrict core motion to the Z- 
axis, 3D computations impose no such constraint and allow 
one to measure the PNS spin in addition to recoil. Recently, 
[Wongwathanarat et al.| ( [2010| ) presented the first 3D excised- 
core calculations of recoil. High PNS velocities were achieved 
for rotating and non-rotating progenitors, providing further 
evidence that PNSs are naturally accelerated during core 
collapse. However, in the case of pulsar spins, [Rantsiou et al.| 
(20111 showed that excising the PNS from the computa- 



tional domain can lead to qualitatively different results in 
the spin rates. As such, future 3D calculations which include 
the PNS should be performed and differences between kicks 
from different progenitor models (rotating and non-rotating) 
should be investigated. 



5 CONCLUSIONS 

We have carried out a suite of axisymmetric simulations 
of the collapse of a massive star's core with the AMR, 
radiation-hydrodynamic code CASTRO. For each calcula- 
tion, we follow the core collapse, PNS formation, ensuing 
neutrino-driven supernova explosion and PNS recoil. By in- 
corporating the effects of neutrino heating and cooling in 
place of more detailed and computationally expensive neu- 
trino transport, we are able to perform multiple calculations 
that simultaneously follow the evolution of the PNS and the 
global explosion for ~1 second and to distances of ~10,000 
km. 

The PNSs in our simulations achieved recoil velocities 
comparable to the those of observed pulsars. After ~1 sec- 
ond of post-bounce evolution, the highest PNS velocity ob- 
tained was 620 kms~^ (model L2.3). After ~0.6 seconds of 
post-bounce evolution, this acceleration was supplied pri- 
marily by the gravitational pull of slow-moving ejecta in 
front of the PNS. This gravitational effect dominates the 
late-time PNS acceleration in all of our calculations. While 
our PNSs have started to decouple from the surrounding 
fluid (see Fig. [3|, the substantial and ongoing gravitational 
acceleration suggests that higher velocities may ultimately 
be achievable. 

Our results suggest that hydrodynamic recoil during an 
asymmetric supernova explosion provides a natural explana- 
tion for pulsar kicks. After the bounce shock stalls, hydro- 
dynamic instabilities deform the shocked material and en- 
sure that the ensuing explosion is asymmetric. Recoil during 
the supernova explosion and gravitational interaction with 
the expanding ejecta subsequently accelerate the PNS to 
high velocities. The results presented in this work are con- 



sistent with the findings ofiNordhaus et al. (2010a I and pre- 



vious axisymmetric simulations that excised the PNS from 
the computational domain ( [Scheck et al.|2004||2006| . Taken 
together, these studies strongly suggest that generic core- 
collapse supernovae can accelerate neutron stars to the high 
velocities observed in the pulsar population. Additionally, 
these studies demonstrate that hydrodynamic processes, and 
not asymmetric neutrino emission, are responsible for this 
acceleration ([Scheck et al.|2006| [Burrows et al.|2007b||Nord-| 
haus et al.|2010a |. In fact, recent simulations of neutron star 



kicks in three dimensions suggest that velocities compara- 
ble to those from axisymmetric calculations are achievable 
( Wongwathanarat et al.|2010 |. 

In this work, we have provided substantial numerical 
support to the hydrodynamic mechanism of pulsar kicks. Re- 
coil due to a neutrino-driven, core-collapse supernova explo- 
sion provides a natural explanation for pulsar kicks without 
appealing to more exotic scenarios. As computational meth- 
ods and resources improve, self-consistent three-dimensional 
studies will enable a full comparison of theoretical models 
to observed distributions of pulsar kicks and spins. 
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